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Abstract. We consider an inertial model of chemotactic aggregation generalizing the Keller-Segel model 
and we study the linear dynamical stability of an infinite and homogeneous distribution of cells (bacteria, 
amoebae, endothelial cells,...) when inertial effects are accounted for. These inertial terms model cells di- 
rectional persistance. We determine the condition of instability and the growth rate of the perturbation 
as a function of the cell density and the wavelength of the perturbation. We discuss the differences be- 
tween overdamped (Keller-Segel) and inertial models. Finally, we show the analogy between the instability 
criterion for biological populations and the Jeans instability criterion in astrophysics. 

PACS. 05.45.-a Nonlinear dynamics and nonlinear dynamical systems 



1 Introduction 

The self-organization of biological cells (bacteria, amoe- 
bae, endothelial cells,...) or even insects (like ants) due 
to the long-range attraction of a chemical (pheromone, 
smell, food,...) produced by the organisms themselves is a 
long-standing problem in physical sciences pP. This pro- 
cess is called chemotaxis. The chemotactic aggregation of 
biological populations is usually studied in terms of the 
Keller-Segel model [2]: 



e^=V-(7^2Vp-i?iVc), 



dc 



-k{c)c + pf{c) + DAc, 



(1) 



(2) 



reproduce the chemotactic aggregation (collapse) of bio- 
logical populations when the attractive drift term Z?iVc 
overcomes the diffusive term D2V p. 

However, recent experiments of in vitro formation of 
blood vessels show that cells randomly spread on a gel ma- 
trix autonomously organize to form a connected vascular 
network that is interpreted as the beginning of a vascula- 
ture [3]. This phenomenon is responsible of angiogenesis, 
a major actor for the growth of tumors. These networks 
cannot be explained by the parabolic model that 
leads to pointwise blow-up. However, they can be recov- 
ered by hyperbolic models that lead to the formation of 
networks patterns that are in good agreement with exper- 
imental results. These models take into account inertial 
effects and they have the form of hydrodynamic equations 



which consists in a drift-diffusion equation ([T]) governing 
the evolution of the density of cells p(r, t) coupled to a dif- 
fusion equation Q involving terms of source and degra- 
dation for the secreted chemical c(r,t). The chemical is 
produced by the organisms (cells) at a rate /(c) and is de- 
graded at a rate k{c). It also diffuses according to Fick's 
law with a diffusion coefficient D. The concentration of 
cells changes as a result of an oriented chemotactic mo- 
tion in a direction of a positive gradient of the chemical 
and a random motion analogous to diffusion. In Eq. ([T|), 
D2{p, c) is the diffusion coefficient of the cells and Di{p, c) 
is a measure of the strength of the influence of the chemi- 
cal gradient on the flow of cells. These coefficients depend 
a priori on the concentration of cells and on the concen- 
tration of the chemical. The Keller-Segel model is able to 



dp ^ 



du 

dc 

di " 



(pu) = 0, 

V)u = --S/p- 
P 



Vc, 



-kc + fp + DcAc. 



(3) 
(4) 
(5) 



The inertial term models cells directional persistance and 
the general density dependent pressure term — Vp(p) can 
take into account the fact that the cells do not interpen- 
etrate. In these models, the particles concentrate on lines 
or filaments |3l4j . These structures share some analogies 
with the formation of ants' networks (due to the attrac- 
tion of a pheromonal substance) and with the large-scale 
structures in the universe that are described by similar 
hydrodynamic (hyperbolic) equations, the Euler-Poisson 
system [5] . The similarities between the networks observed 
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in astrophysics (see Figs 10-11 of [6]) and biology (see Figs 
1-2 of ^) are striking. 

In order to make the connection between the parabolic 
model (II])-(l2|) and the hyperbolic model ([S])-©, we con- 
sider a model of the form 



dp „ 



(pu) = 0, 



(6) 



2 Instability criterion for biological 
populations 

2.1 The dispersion relation 

We consider an infinite and homogeneous stationary so- 
lution of Eqs. (H])-® with u = 0, p = Cst. and c = Cst. 
such that 



^ + (u.V)u 



= -DaVp + L>iVc - ^pu, (7) 



= /(c)p. 



(10) 



p/(c) + DAc, 



(8) 



Linearizing the equations around this stationary solution, 
we get 



including a friction force — ^pu. This type of damped hy- 
drodynamic equations was introduced in [718] at a general 
level. This inertial model takes into account the fact that 
the particles do not respond immediately to the chemo- 
tactic drift but that they have the tendency to continue in 
a given direction on their own. However, after a relaxation 
time of order their velocity will be aligned with the 
chemotactic gradient. This is modeled by an effective fric- 
tion force in Eq. ([7]) where the friction coefficient ~ 
is interpreted as the inverse of the relaxation time. This 
term can also represent a physical friction of the organ- 
isms against a fixed matrigel. In the strong friction limit 
^ — > -|-oo, or for large times t ^ C"^? ons can formally 
neglect the inertial term in Eq. ((T]) and obtain j8|9|10j : 



(11) 



p^ = -D2V5p + DiV5c-ip5u, (12) 
ot 



= (/'(c)p - k)5c + f{c)Sp + DASc, (13) 

ot 

where we have set k ~ k{c) + ck'(c). Eliminating the ve- 
locity between Eqs. pT]) and IT^ . we obtain 



d'^Sp 86 p 



dt^ ^ dt 



^—^ = D2ASP- DiASc. 



(14) 



pu 



-^{D2Vp 



DiVc) + 0(r 



(9) 



Looking for solutions of the form dp ~ 5 pe"*' e^"^'^ and 5c ■ 
Sce'^^e^^ '^ , we get 



Substituting this relation in Eq. we recover the parabolic 
Keller-Segel model Therefore, the Keller-Segel model 

can be viewed as an overdamped limit of a hydrodynamic 
model involving a friction force. Alternatively, neglecting 
the friction force ^ = 0, we recover the hydrodynamic 
model introduced in [3]. 

In this paper, we study the linear dynamical stability of 
an infinite and homogeneous distribution of cells with re- 
spect to the inertial model (H])-®. We determine the con- 
dition of instability and the growth rate of the perturba- 
tion, and discuss the differences with the results obtained 
with the Keller-Segel model ([I])-©. We also discuss some 
analogies with the dynamical stability of self-gravitating 
systems. Indeed, there are many analogies between the 
chemotactic aggregation of biological populations and the 
dynamics of self-gravitating Brownian particles [9| . In par- 
ticular, the Keller-Segel model (P)-© is similar to the 
Smoluchowski-Poisson system [11] and the hydrodynamic 
equations ([6])- ((HI) are similar to the damped Euler equa- 
tions of Brownian particles [8 lOj . The main difference 
between biological systems and self-gravitating Brownian 
particles is that the Poisson equation in the gravitational 
problem is replaced by a more general field equation ([H]) 
taking into account the specificities of the biological prob- 
lem. Owing to this analogy, we shall discuss the relation 
between the instability criterion of biological populations 
and the Jeans instability criterion T2] in astrophysics. 



{F - ea)Sc + f{c)Sp = 0, (15) 



DiqHc-iD^q^ +aia + O)Sp = 0, (16) 

where F = /'(c)p — k — q^D. These equations have non- 
trivial solutions only if the determinant of the system is 
equal to zero yielding the dispersion relation 

ea^ + (e^ - F)a^ - {F^ - eq^D^)^ 

-q^{f{c)Di+D2F)^Q. (17) 

The condition of marginal stability {a = 0) corresponds 
to f{c)Di + D2F = and the condition of instability is 



f{c)Di+D2F > 0. 



(18) 



We note that the instability criterion does not depend 
on the value of e and ^. Let us consider in detail some 
particular cases. 



2.2 The case e = 

When the chemical has a large diffusivity, the temporal 
term in Eq. ^ can be neglected fT3\. Thus, we formally 
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consider e = 0. In that case, the dispersion relation ([T7] 
reduces to 



2 , e , 2 / n . /W-^i 



F 



The discriminant is 



A{q) =e- V D2 + 



f{c)Di 
Fiq) 



and the two roots are 



(19) 



(20) 



(21) 



If Reijj) < the perturbation decays exponentially rapidly 
and if Re{(j) > the perturbation grows exponentially 
rapidly. In that case, the system is unstable and Reicr) is 
the growth rate of the perturbation, li D2+ f{c)Di/ F < 0, 
then Z\ > so that cr+ > (unstable). Alternatively, if 
D2 + f{c)Di/F > 0, then A < C^. Either Z\ < and 
i?^(cr) = -C/2 < (stable) or < Z\ < ^nd cr± < 
(stable). Therefore, the system is unstable if 



D2<^ 



k- f'{c)p + Dq^ 



(22) 



and stable otherwise. To determine the range of unstable 
wavelengths, we must consider different cases: 



2.2.1 If fc-/'(c)p> 0: 

In that case, a necessary condition of instability is that 



D2<=^ 



k - J'{c)p 



^ {D2)c 



(23) 



If this condition is fulfilled, the unstable wavenumbers are 
determined by 



f{c)Di 
D2 



+ f'{c)p~k 



_ 2 

— ^max ■ 



(24) 




Fig. 1. Graphical construction determining the range of un- 
stable wavenumbers. We have taken ^ = D2 = D = 1, fDi = 2 
and k — f' {c)p = 1 (solid line) or k — f' {c)p — (dashed line). 



0.15 - 



o 0.1 - 




0.05 - 



Fig. 2. Evolution of the growth rate of the perturbation as a 
function of the wavenumber for different values of e. We have 
taken ^ = D2 = D = 1, fDi = 2 and - f'{c)p = 1. 



The growth rate of the perturbation with wavenumber q is jj-^ particular for ^ = we have 
a+ = ^{—£_ + ^/A{q)). Therefore, the most unstable mode 
(7* is the one which maximizes A{q). It is given by 



Dql 



f{c)D,(k-f{c)p) 



Do 



1/2 



+ f{c)p-k. (25) 



2(7* = 



4/(c)i?i 



D 



1 - 



iD2{k-r{c)p) 



(27) 



The largest growth rate cr, is then determined by 



2a, = 



4.f(c)ZJi 
D 



1 - 



lD2ik-f'ic)p) 



f{c)Di 



(26) 



The range of unstable wavelengths is determined graphi- 
cally in Fig. [T] and the evolution of the growth rate of the 
perturbation as a function of the wavenumber is plotted in 
Fig. [2] (we have also consider the case e ^ in this Figure 
by solving Eq. pT| which is a second degree equation in 
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2.2.2 If fc-/'(c)p = 0: 

In that case, {D2)crit - 
are determined by 



-oo. The unstable wavenumbers 



(28) 



The most unstable mode is = and the largest growth 
rate cr* is given by 



2cr. 



4/(c)gi 
D 



(29) 



2.2.3 If fc-/'(c)p< 0: 

In that case, the system is unstable for 



f'{c)p-k 
D 



2 1 



f{c)D 



Do 



l + f'{c)p-k 



The growth rate diverges when 



(30) 



(31) 



corresponding to F(q) — 0. Close to the critical wavenum- 
ber go J we have 



2D 



1 



(q ^ q+). (32) 



This expression is valid for ^ finite and e = 0. It can be 
directly obtain from Eq. by using F ^ ~2DqQ{q — 
9o) for q qo- Thus, when the temporal term is 
neglected in Eq. ([5]), i.e. e = 0, a critical behaviour occurs. 
This critical behaviour is regularized for e ^ 0. Indeed, 
taking q = go, i-e. F ^ 0, in Eq. (|17[) . we obtain 



ea^ + eia^ + eq^D^a ~ qlf{c)D^ - 0. 
Taking the limit e — )■ 0, we find that 

1/3 



^(go) 



4fic)D, 



(33) 



(34) 



which is finite for e > but diverges like when e ^ 0. 
For e ^ and q ~^ qo, the dispersion relation can be 
simplified in 



ea^ + 2Dqa{q - qa)a^ - q^f{c)D, - 0. 



(35) 



For e = we recover Eq. p2p and for q = qo we recover 
Eq. p4p . For q ^ qo, we can easily express g as a function 
of a according to 



qlf{c)D, 



ecr 



(36) 



2L>goCT2 

On the other hand, for g = 0, Eq. p7|) reduces to 

e^' + (ee - /'(c)p + k)cj^ - af{c)p ~ k)a - 0. (37) 











Unstable 










Pmax 







Fig. 3. Graphical construction determining the range of un- 
stable wavenumbers. We have taken ^ = D2 = D = 1, fDi — 2 
and k — f'{c)p = —1. 



The positive root of this equation is 
a(0) = = 



(38) 



which is finite for e > but diverges like when e — > 0. 

The range of unstable wavelengths is determined graph- 
ically in Fig. [3] and the evolution of the growth rate of the 
perturbation as a function of the wavenumber is plotted in 
Fig. m (we have also consider the case e 7^ in this Figure 
by solving Eq. ([T7]))- We see that the case e = is very 
special. For e = 0, the range of wavenumbers q < qo seems 
to be stable according to Eq. ((30)) because the two roots 
of Eq. (fT9|) are negative. However, for any finite value of 
e, a third root appears. This root is positive and tends to 
infinity when e ^ (see Eqs. (p4)) and p8|) ). Therefore, 
this unstable branch is rejected to infinity when e —^ 0. 
This implies that the region g < go is in fact extremely 
unstable for e = 0^. In particular, for e > 0, the most 
unstable mode is g, = and the largest growth rate is 
given by Eq. 



-I-CXD 



2.3 The case ^ 



In the overdamped limit, the hydrodynamical equations 
dH)-® return the Keller-Segel model H])-®. Let us con- 
sider the stability analysis in that case for comparison with 
the inertial case. The dispersion relation now reads 



e^a^ - (FC - eq^D2)a - q\f{c)Di 
and the two roots are 



D2F) ^ 0, (39) 



CT± 



f C - e^2g^ ± ^Ajq) 
2,i 



(40) 



where 



^(9) = [n + fq^D2f + Ae^q^f{c)Di > 0. (41) 
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Fig. 4. Evolution of the growth rate of the perturbation as a 
function of the wavenumber. We have taken ^ = D2 = D — 1, 
fDi =2 andl-/'(c)p= -1. 



6 



<T 4 



Unstable 



f(c)D,/D, 



Fig. 5. Graphical construction determining the range of un- 
stable wavenumbers. We have taken D2 = D — 1, fD\ — 2 
and k — f'{c)p = 1 (solid line) and fe — /'(c)p = — 1 (dashed 
line) . 



Writing the solution in the form 
-b±Vl^ 



Aac 



2a 



(42) 



it is easy to see that the system is stable if (i) F < eq^D2/£, 
{b > 0) and if (ii) F < -f{c)Di/D2 (c > 0). It is unstable 
otherwise. Since (ii) implies (i), the system is stable if F < 
—f{c)Di/D2 and unstable otherwise. Thus, the system is 
unstable if 



f'{c)p + Dq^< 



f{c)Dl 
D2 



(43) 



and stable otherwise. The range of unstable wavelengths 
is determined graphically in Fig. El 
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2.3.1 If fc-/'(c)/9> 0: 
In that case, the system is unstable for 

/(c)i^i 



D2 < 



k- f'{c)p + Dq^' 



(44) 



and stable otherwise. This is the same criterion as for the 
inertial model. A necessary condition of instability is that 



k-r{c)p ^ 



(45) 



If this condition is fulfilled, the unstable wavenumbers are 
such that 



7 1 



f{c)D 



Do 



l + f'{c)p-k 



= 9m, 



(46) 



Fig. 6. Evolution of the growth rate of the perturbation as 
a function of the wavenumber for the Keller-Segel model. We 
have taken D2 = D = 1, fDi = 2 and I - /'(c)p = 1. 



We now determine the value of the optimal (most un- 
stable) wavenumber and the corresponding growth rate 
CT* (see Fig. Eq. is a second order equation in cr 
whose solutions are given by Eq. (I40p . We can maximize 
<y+{q) to obtain q^ and cr*. However, it appears simpler 
to proceed differently. Eq. ((39)l can also be viewed as a 
second order equation in a; = of the form 



with 



Ax^ + B{a)x + C{a) = 0, 



A^DD2, 



(47) 
(48) 



These results are unchanged with respect to the inertial 
case. 



B{a) - {D^ + eD2)a - f{c)D^ - D2{f{c)p - k), (49) 
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C{a)^e^a'-af{c)p-ky>0. 



(50) 



There will be two roots xi and X2 provided that B(a) < 
and Z\(cr) ^ - A AC > 0. This last condition can be 
written 



A{a) 



ba + c>0, 



with 



(51) 
(52) 



b = ^2[D2{fic)Di + D2if'ic)p - k))e 
+DCfic)Di - DD2af{c)P - k)] < 0, (53) 

c^[fic)D,+D2inc)p-k)f. (54) 

The discriminant S = ~ Aac of Eq. (|5ip is given by 

S = 16fic)^DDiD2[{fic)D, 
+D2{f'{c)p - k))e - Danc)p - k)]. (55) 

The condition A{(7) > to have two roots xi and X2 is 
equivalent to cr < cr* with 



-b^VS 
2a ' 



(56) 



(Note that the possibility a > {—b + ^f5)/2a must be 
rejected since it does not satisfy the requirement B{a) < 
0). For (T = cr, , the two roots xi = X2 = x^ coincide 
(A = 0) so that cr* is the maximum growth rate. It is 
reached for an optimal wavenumber — —B/(2A), i.e. 



-B{a,) 
2A ' 



(57) 




Fig. 7. Evolution of the growth rate of the perturbation as 
a function of the wavenumber for the Keller-Segel model. We 
have taken D2 = D = 1, fDi = 2 and - f'{c)p = -1. 



2.4.1 \f k - f{c)p>0: 



This is a particular case of Sec. 12.2.11 corresponding to 
^ +00. The expression of the largest growth rate is 
given by 



1 



iD2{k-r{c)p) 



The other results are unchanged. 



(61) 



2.3.2 If fc - f'{c)p < 0: 2.4.2 H k - f'{c)p < 

In that case the system is unstable for the wavenumbers 
1 



+ / (c)p - k 

U2 



2 



Imax- (58) 



The system is unstable for the wavenumbers determined 
by Eq. ([50]) . For , corresponding to F ~ —2Dqo{q— 

Qo) 0, the growth rate diverges like 



The growth rate of the perturbation as a function of the qof{c)Di 

wavenumber is plotted in Fig. [71 As discussed in Sec. 12. 2. 3] ~ 2D{q — qo) ' 
the case e = is special and will be considered specifically 

in the next section. This divergence is regularized if e 7^ 0. Taking q = qo, i.e. 

F = in Eq. we get 

2.4 The case f ^ +00 and e = g^^2 ^ ^q^j^^^ _ qlf{c)Di = 0. (63) 

If we neglect the temporal term (e — 0) in the Keller-Segel e > we obtain 

model (^ — !■ +00), we obtain the dispersion relation 

F^a + q\f{c)D, + D2F) = 0, (59) ^^^^^ ^ f <llf{c)DA ^^^^ 



so that a is explicitly given by 



2 / f{c)Di \ which is finite for e > but diverges like e when e — > 0. 

= 9 yjjq2^k - f'{c)p ■ ^^^^ For e ^ and q q+ , Eq. §^ can be simplified in 

The instability criterion is given by Eq. e^cr^ + 2Dqo({q - qo)<7 ~ qlf{c)Di = 0. (65) 
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For e = we recover Eq. (|62p and for q — qo we recover 
Eq. dMl)- The solution of Eq. §Bi is 



ecr = -Dqa{q - qo) + W D^q^{q - qof + 



On the other hand, for (7 = 0, Eq. (1551) leads to 

.(0) = l^^^, 



(66) 



(67) 



which is finite for e > but diverges hke e ^ when e — > 0. 

Equations ([221)- (EH) and Eqs. ([121)- dMl) differ because 
both (7 and ^ tend to infinity, so that the expression de- 
pend on how the Hmits are taken. The general case can be 
treated as follows. For e = 0, taking g — > in Eq. ^ 
we get 



<y{'y + e) 



go/(c)Ji 
2D{q - qo) 



(68) 



On the other hand, for e 7^ 0, taking q — qo (i-C. F — 0) 
in Eq. (fTT)) . we get 



90/(0)^1 



Finally, for e — )■ and q ^ qo, we have 

eCT3 + (e^ + 2Z?(7o((?-(?o))'T2 
+2Dqo{q - qo)((J - qlf{c)Di - 0. 



which reproduces the correct behaviours 



(69) 



(70) 
and (pg)) . 



3 Analogy with the Jeans problem in 
astrophysics 

3.1 The damped Euler equations 

Let us consider a particular case of Eqs. (H])-® corre- 
sponding to D2{p, c) — p'{p) and Di{p, c) — pS'{c) where 
p and S are arbitrary functions. In that case, the hydro- 
dynamical equations take the form 



dp 
dt 



V • (pu) = 0, 



(71) 



^-l-(u.V)u = — Vp + V5(c)-^u, (72) 
at p 



dc 

^7i7 = -k{c)c + pf{c) + DAc. 
at 



(73) 



For ^ +00, we can neglect the inertia of the particles 
so that pu ~ — i(Vp — /9VS'(c)). Substituting this relation 



in Eq. (j7ip . we obtain a special case of the Keller-Segel 
model 



dp 
dt 



V-[x {Vp-pS'{c)\/c)] 



(74) 



where we have set x = l/'C- Equations ((7T|) - ((72)) can be 
viewed as fluid equations appropriate to the chemotactic 
problem. Equation (j7ip is an equation of continuity and 
Eq. ((72|) is similar to the Euler equation where p plays the 
role of a pressure and the chemotactic attraction plays the 
role of a force. Since p = p{p), these equations describe a 
barotropic fluid. The main novelty of these equations with 
respect to usual hydrodynamical equations is the presence 
of a friction force which allows to make a connection be- 
tween hyperbolic (^ = 0) and parabolic (^ —t +00) mod- 
els. 

This hydrodynamic model including a friction force is 
similar to the damped barotropic Eulcr-Poisson system 
which describes a gas of self-gravitating Brownian parti- 
cles |8ll0j or the violent relaxation of coUisionless stellar 
systems on the coarse-grained scale in astrophysics [7] . In 
that analogy, the concentration of the chemical c plays the 
role of the gravitational potential <P. The main difference 
between the two models is that the Poisson equation for 
self-gravitating systems is replaced by a more general field 
equation (j73p for bacterial populations. To emphasize the 
connection with astrophysical problems, let us consider a 
particular case of Eqs. ((7T|) - ([73|) where e = 0, S{c) — c and 
k and / are constant. Then, introducing notations similar 
to those used in astrophysics (noting c = — k/D 
f/D = SdG), we can rewrite Eqs. ^-1^ in the form 



i-2 



| + V.(pu).0, 



du , , 1 

— + (u • V)u = --Vp - V<? - ^u, 



kl<P 



SdG{p - p). 



(75) 



(76) 



(77) 



When fco = P = 0, these equations are isomorphic to the 
damped Eulcr-Poisson system describing self-gravitating 
Brownian particles [8|10j . In the strong friction limit, we 
get 



dt 



i (Vp + pW-P) 



(78) 



which can be interpreted as a generalized Smoluchowski 
equation. Thus, for ^ — > -l-oo, we obtain the generalized 
Smoluchowski-Poisson system describing self-gravitating 
Brownian particles in an overdamped limit [8]. Alterna- 
tively, for ^ = we recover the barotropic Euler-Poisson 
system that has been studied at length in astrophysics to 
determine the period of stellar pulsations [Uj and the for- 
mation of large-scale structures in cosmology [S]. There- 
fore, the chemotactic model ((75|) - (|77p is similar to astro- 
physical models with additional terms. In the astrophysi- 
cal context, the case fco 7^ would correspond to a shield- 
ing of the gravitational interaction on a typical length k^^ . 
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This Yukawa shielding appears in theories where the gravi- 
ton has a mass but in that case fco is very small which does 
not need to be the case in the biological problem. 

We now consider the linear dynamical stability of an 
infinite and homogeneous solution of Eqs. ([75)1 - ([T7)l . By 
mapping the equations (|75|) - ((77|) onto a generahzed as- 
trophysical model, we shall see that the instability crite- 
ria obtained in Sec. are connected to (and extend) the 
Jeans instability criterion of astrophysics [T2]. To avoid 
the Jeans swindle (15^ when ko = 0, we have introduced 
a "neutralizing background" p in Eq. (|77p. In fact, in the 
biological problem, this term appears naturally when we 
consider the limit of large diffusivity of the chemical (see 
[13]); therefore, there is no "Jeans swindle" in the biolog- 
ical problem based on Eq. A similar term p appears 
in cosmology when we take into account the expansion of 
the universe and work in the comoving frame [5]. 



3.2 The dispersion relation 

We consider an infinite and homogeneous stationary solu- 
tion of Eqs. (Hgi-dTTI) with u = 0, p = Cst. and <P = Cst. 
such that 



(79) 



For a pure Newtonian interaction with fcg = 0, we have 
p = p. Linearizing the equations around this stationary 
solution, we get 



d6p 



pV • (Su = 0, 



dSu 



-c^Vf^p - pV5<P - ^p5u, 



A6<P - klS<P = SdGSp, 



(80) 

(81) 
(82) 



where we have introduced the velocity of sound cj. = p'{p). 
Eliminating the velocity between Eqs. and (|5T|) . we 
obtain 



d'^Sp dSp 



= ci^Sp + pAS^. 



(83) 



dt^ ^ dt 

Looking for solutions of the form 5p ^ (Jpe'^'' ''""*', we get 
(-^2 - i^u + clk^)5p = -pk^SS, (84) 



SdG ^ 
u2 TT^—.^P- 



(85) 



k'^ + k^ — iuj 

From these equations, we obtain the dispersion relation 



uj{uj + iO ^ cik 



2, 2 _ SdGpk'^ 
k^ + kf 



(86) 



In the case ^ = and feg = 0, we recover the usual Jeans 
dispersion relation [T?| : 



u;2 = clk^ ^ SdG p. 



(87) 



3.3 Instability criterion 

If we set a — —iuj, the dispersion relation becomes 

SdGp 



a^+^tJ + k^ 



0. 



The two roots are 



with 



Accordingly, the system is unstable if 

SdGp 



k^ + kl' 



(89) 



(90) 



(91) 



and stable otherwise. A necessary condition of instability 
is 



ct < 



SdGp 



(92) 



If this condition is fulfilled the unstable wavelengths are 
such that 



, ^ / SdGp 2 _ I 
- \ r.2 - f^r. 



(93) 



The wavelength which has the largest growth rate is given 

by 



,_( SdGpkl V'' 
and the corresponding growth rate is given by 



(94) 



2(7 * 



e + 45rfGp 1 



SdGp 



(95) 



The instability criterion (|9ip is equivalent to the insta- 
bility criterion ((22|) of Sec. 12.21 for the particular case of 
the chemotactic model considered in Sec. 13.11 The parallel 
with astrophysics is interesting to develop in order to give 
a more physical interpretation to Eq. (|22p. All the other 
formulae can be interpreted accordingly. In particular, we 
note that the coefficient D2 plays the role of a velocity of 
sound cl. 



3.4 Particular cases 

Let us consider particular cases of the foregoing expres- 
sions: 
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T<Tc 


Unstable \ 













1 2 3 4 5 

k/k„ 

Fig. 8. Graphical construction determining the range of un- 
stable wavenumbers. 




• For = 0, one has A:„ 



-co, fc, — +00 and 



Fig. 9. Maximum wavenumber kmax{T) and most unstable 
wavenumber kt{T) as a function of the temperature T. The 
line kmax{T) determines the separation between stable and 
unstable states. 



(96) 



• For fco = (Newtonian potential), one has {0^)0^1 



-00, k„ 



and 



1 /"^ 

[SdGp/cj.) " = kj (Jeans length), k^ = 



For ^ = 0, one has 



For ^ +00, one has 



SdGp 
cr* = — - — I 1 



SdGp 



(97) 



(98) 



(99) 



the growth rate of the perturbation can be written 



'1 - 



4SdGpk^ (T 



C2 Ic^ \ T 



1 



1 + (fc/fco) 



The condition of instability reads 

T 1 
— < — 



(fc/fco)2 



(101) 



(102) 



and a necessary condition of instability is T < Tc- For 
T < Tc the unstable wavenumbers (see Fig. [5]) are such 
that k < kmax{T) with 



ko 



T 



(103) 



The wavenumber with the largest growth rate is given by 



k.{T) 



1/2 



1/2 



(104) 



3.5 Isothermal gas 

For an isothermal gas with an equation of state p — pT 
(for simplicity, we have noted T instead of kBT/m), the 
velocity of sound is equal to the square root of the tem- 
perature: = T. It is relevant to re-express the previous 
relations as follows. Introducing the critical temperature 



SdGp 



fr2 



(100) 



and the largest growth rate by 



T 



-1 + 



'^SdGp 

e 



1/2 



(105) 



The maximum wavenumber k„iax (T) and the most unsta- 
ble wavenumber k^{T) are plotted as a function of the 
temperature in Fig. [S) In Fig. llOl we represent the growth 
rate a{k) as a function of the wavenumber. Finally, in Fig. 
[Tl] we plot the largest growth rate cr»(T) as a function 
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0.03 



0.02 



0.01 



T/Tc=0.5 






^^^^^^^ a.(T),k.(T) 




k™<(T)\ 



0.2 



0.4 0.6 

k/k„ 



Fig. 10. Growth rate of the perturbation as a function of the 
wavenumber. We have taken ASdCp/^"^ = 1 and T/Tc = 0.5. 




Fig. 11. Dependence of the largest growth rate o-,{T) with 
the temperature. We have taken 4SdGp/^^ = 1. 



of the temperature. The maximum value of the largest 
growth rate cr, (T) is obtained for T = and is given by 



2(cr*)f 



1 + 



iSdGp 



(106) 



For fco = (Newtonian interaction), then Tc = +oo 
and the growth rate of the perturbation can be expressed 
as 



2cr 

T 



(107) 



The condition of instablity is 



SdGjA 
T 



1/2 



(108) 



0.4 



0.3 



0.2 



0.1 



0.2 



0.4 0.6 

k/k„a,(T) 



Fig. 12. Growth rate cr(fc) as a function of the wavenumber 
when ko = 0. We have taken iSdGp/^^ = 1. In terms of the 
variable k / kmax(T) , the curve is independent on the tempera- 
ture. 



where k^ax is the equivalent of the Jeans wavenumber. 
The maximum growth rate is obtained for /c* = and its 
value is given by 




(109) 



independently of the temperature. The growth rate cr(fc) 
is represented as a function of the wavenumber in Fig. 1121 



4 Conclusion 

In this paper, we have studied the linear dynamical sta- 
bility of an infinite and homogeneous distribution of bio- 
logical cells whose density distribution evolves under the 
process of chemotaxis. We have modeled their evolution by 
hydrodynamical equations including a friction force [8l9j . 
This inertial model takes into account the fact that the 
cells do not respond immediately to the chemotactic drift 
but that there is a relaxation time for their velocity 
to get aligned with the chemotactic gradient. The usual 
Keller-Segel model [2] is recovered in the strong friction 
limit ^ — » -|-oo (or for large times t ^ Alternatively, 
for ^ = 0, we recover the inertial model of [3]. We have 
shown that these equations were similar to those describ- 
ing self-gravitating Brownian particles and that the dy- 
namical stability of biological populations was related to 
the Jeans problem in astrophysics. These results extend 
the analogies between biology and astrophysics investi- 
gated in [3. 

The mathematical model (H])-® considered in this pa- 
per can have applications in biology. Depending on the 
value of the parameters, it can decribe different sorts of 
systems. The Keller-Segel model ([I])-© obtained in the 
overdamped limit ^ — > c» in which inertial terms can be 
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neglected is appropriate to describe experiments on bacte- 
ria like Escherichia Coli and slime mold amoebae like Dic- 
tyostelium discoideum 2^ . These systems exhibit pointwise 
concentrations as a result of chemotactic collapse. On the 
other hand, the hydrodynamic model ©-(HI) was shown 
to generate a vascular network starting from randomly 
seeded endothelial cells [3 4'. This can account for exper- 
iments of in vitro formation of blood vessels where cells 
randomly spread on a gel matrix autonomously organize 
to form a connected network that can be interpreted as 
the initiation of angiogenesis. This is also similar to the 
formation of capillary blood vessels in living beings dur- 
ing embriogenesis [TO . The authors of [3] evidence a per- 
colative transition as a function of the concentration of 
cells. Above a critical density, the system forms a contin- 
uous multi-cellular network which can be described by a 
collection of nodes connected by chords. For even higher 
concentrations a "swiss cheese" pattern is observed. Such 
structures can be obtained only if inertial terms are ac- 
counted for. Another hydrodynamic model of bacterial 
colonies taking into account interial terms has been pro- 
posed by Lega & Passot [T7] to describe the evolution of 
bacterial colonies growing on soft agar plates. This model 
consists in advection-reaction-diffusion equations for the 
concentrations of nutrients, water, and bacteria, coupled 
to a single hydrodynamic equation for the velocity field 
of the bacteria-water mixture. This model is able to re- 
produce the usual colony shapes together with nontriv- 
ial dynamics inside the colony such as vortices and jets 
recently observed in wet colonies of Bacillus subtilis [TH]. 
This can be linked to a process of inverse cascade of energy 
as in two-dimensional hydrodynamic turbulence. Lega & 
Passot [T7] show that the large-scale Reynolds numbers 
can be relatively high so that inertial effects have to be 
taken into account to adequately model the experiments 
of [18] . It is shown also that viscosity is important in this 
model. Although the hydrodynamic equations of [17] are 
different from Eqs. (I6|)-([8|), their model displays a mech- 
anism for collective motion towards fresh nutrients which 
is similar to classical chemotaxis. In particular, a chemo- 
tacticlike behaviour and a connection to the Keller-Segel 
model (HI-© is obtained for short times. 

In this paper, we have considered solutions of Eqs. 
(©-([H) near an infinite and homogeneous distribution and 
we have investigated the time dependence of these solu- 
tions in the linear regime Q- When the criterion (fT5)) is 
fulfilled, the appearance of a spontaneous perturbation 
can lead to an instability. The perturbation grows until 
the system can no longer be described by equilibrium or 
near-equilibrium equations. In that case, we must account 
for the full nonlinearities encapsulated in Eqs. (©-([H]). Of 



^ The linear dynamical stability of inhomogeneous distribu- 
tions of bacteria has also been studied in [lll20j for over- 
damped models and in pIF for inertial models, when the equa- 
tion for the concentration of the chemical takes the form of 
a Poisson equation (|lll|l like in gravity. In these studies, the 
distribution of particles is self-confined [10] or confined in a 
finite domain (box) [11120] . In biology, the box can represent 
a droplet or the container itself. 



course, the nonlinear regime of instability is the most rele- 
vant for biological applications. This nonlinear regime has 
been investigated in detail for a reduced version of the 
Keller-Segel model [13] : 



Ac = -\p. 



(pVc), 



(110) 



(111) 



In that case, the concentration of the chemical is related 
to the concentration of the bacteria by a Poisson equa- 
tion. These equations are isomorphic to the Smoluchowski- 
Poisson system 



A$ = SdGp. 



(112) 



(113) 



describing self-gravitating Brownian particles '11]. In di- 
mension d > 2, they exhibit blow-up solutions leading ulti- 
mately to the formation of Dirac peaks. This corresponds 
to a chemotactic collapse in biology or to an isothermal 
collapse (in the canonical ensemble) in gravity. There is a 
vast literature on the theoretical study of these equations 
both in applied mathematics (see the review by Horstmann 
[I9] ) and in physics [11|20|21|22|T0] . Generalized chemo- 
tactic models and generalized gravitational models have 
also been studied, like in [23 to account for anomalous dif- 
fusion or like in [1^ to account for inertial effects. On the 
other hand, bifurcations between "stripes" and "spots" 
have been found when the degradation of the secreted 
chemical is taken into account so that the equilibrium 
structures of the bacterial colonies are similar to "domain 
walls" in phase ordering kinetics [24]. The linear insta- 
bility regime that we have considered in this paper initi- 
ates the nonlinear regime where interesting and non-trivial 
structures form, accounting for the morphogenesis of bac- 
terial populations. In the linear instability analysis, the 
general form of perturbation is a superposition of sinu- 
soidal waves. Each single wave corresponds to a "streak" 
with relatively high density. However, other patterns like 
regularly spaced "clouds" can be obtained by a proper su- 
perposition of "streaks" (see Appendix A of |2]). These 
"clouds" will be presumably selected by nonlinear effects 
and each of them can initiate a local collapse leading to 
pointwise blow-up |19|11|20] . Indeed, these clouds have 
the radial symmetry that is assumed at the start in most 
studies of chemotactic collapse. This will lead to a set of 
N singular structures. These compact structures interact 
with each other and lead to a coarsening process where the 
number of clusters decays in time as they collapse to each 
other. This process may share some analogies with the 
aggregation of vortices in two-dimensional decaying tur- 
bulence [25]. Therefore, the connection between the linear 
regime investigated in this paper and the nonlinear regime 
investigated in [19|11|20] is relatively clear. 
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